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Abstract. In this paper we present mathematical approaches to understand a sym- 
metry break and formation of spatially heterogenous structures during development. 
We focus on the models given by reaction-diffusion equations and approach the ques- 
tion of possible mechanisms of development of spatially heterogeneous structures. We 
discuss two mechanisms of pattern formation: diffusion-driven instability (Turing in- 
stability) and a hysteresis-driven mechanism, and demonstrate their possibilities and 
constraints in explaining different aspects of structure formation in cell systems. De- 
pending on the type of nonlinearities, we show the existence of Turing patterns, the 
maxima of which may be of the spike or plateau type, and the existence of transition 
layer stationary solutions. These concepts are discussed on example of morphogenesis 
of the fresh water polyp Hydra, which is a model organism in developmental biology. 
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8.1.1 Introduction 


Spatial and spatio-temporal structures occur widely in physics, chemistry and biology. 
In many cases, they seem to be generated spontaneously. Understanding the principles 
of development and design in biology is among the crucial issues not only in devel- 
opmental biology but also in the field of regenerative medicine. In order to develop 
methods for intelligent engineering of functional tissues, the main principles of devel- 
opment and design have to be understood. 

Recent advances made in genetic and molecular biology have led to detailed de- 
scriptions of a number of events in embryological development. Although genes con- 
trol pattern formation, genetics alone is insufficient to understand which physio-chemi- 
cal interactions of embryonic material produce the complex spatio-temporal signaling 
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cues which ultimately determine the cell’s fate. Since the establishment of symmetry 
breaking by cell polarity in developing tissues is determined by quantitative integra- 
tion of multiple signals in a highly dynamical and self-organized process, it can be 
hardly understood using conventional molecular biology methods alone. The role of 
mathematical modeling is to verify which processes are sufficient to produce the pat- 
terning. Model mechanisms can suggest to the embryologist possible scenarios as to 
how, and sometimes when, a pattern is laid down and how the embryonic form might 
be created. Modeling also allows to make experimentally testable predictions and may 
provide alternative explanations for the observed phenomena. 

In other areas of biology, such as neurophysiology or ecology, mathematical mod- 
eling has led to many discoveries and insights through a process of synthesis and in- 
tegration of experimental data, see [37] and references therein. Also in developmental 
biologymathematical models in developmental biology many different morphologies 
have been the subject of mathematical modeling. Some of the biological systems have 
attained the status of a paradigm in theoretical work [3,7]. 

One such example, which shows how the study of model mechanisms can sug- 
gest real scenarios for the process of pattern formation, is limb development [37]. A 
mechanochemical model describes the diffusion, haptotaxis and advection of mes- 
enchymal cells which evolve in a developing limb bud and which eventually become 
cartilage. The other developmental process for studying pattern formation and differ- 
ent aspects of embryogenesis is the segmentation of the insect embryo [3, 12]. The 
models based on chemotaxis and the response of cells to gradients in the chemoattrac- 
tant were applied to study the life cycle of the slime mould Dictyostelium discoideum 
and emergence of concentrated patterns of cell density [11, 13,42]. More recently, 
models of morphogenesis have been applied to understand the growth of tumors. They 
involve a wide range of biological phenomena such as cell-adhesion and cell traction, 
angiogenesis, pattern formation in cancer and macrophage dynamics [5, 37]. 

All these models, although based on different biological hypotheses, have many 
common mathematical features and are mostly based on a few views of pattern gen- 
eration. One is the chemical prepattern approach involving hypothetical chemicals 
(morphogens) which diffuse and react in sucha way that spatial heterogeneous patterns 
can evolve from the uniform steady states. Coupling diffusion process of signaling 
molecules with nonlinear dynamics of intracellular processes and cellular growth and 
transformation leads to receptor-based models, which differ from the usual reaction- 
diffusion systems. Next, the mechanochemical approach takes into account mechani- 
cal forces and properties of cells and tissues. Another class of models rely on taxis, e.g., 
chemotaxis or haptotaxis, and the response of cells to gradients in the concentration 
of signaling molecules in the environment [1, 30]. 

Different models are able to produce similar patterns. The question is how to dis- 
tinguish between them so as to determine which may be the relevant mechanisms. 
Of course the first necessary condition is that the model must produce observed pat- 
terns. But then it is important to design new experiments, which could allow for model 
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validation and therefore would also lead to the verification of model hypotheses and 
biological theories. 

Whereas different developmental processes are involved in different organisms, it 
is striking how conserved the processes are across different taxonomic units such as 
the phyla. Also, the same processes are involved in various diseases, in particular in 
cancer. Molecules found to be oncogenic factors also play an important role in devel- 
opmental processes. This unity and conservation of basic processes implies that their 
mathematical models can have an impact across the spectrum of normal and patholog- 
ical development. 

In this paper we present different approaches to model development and regenera- 
tion of the fresh-water polyp Hydra. We focus on the framework of continuous mod- 
els given by partial differential equations. In particular, we employ reaction-diffusion 
equations and discuss their performance in the context of key experiments. One of the 
objectives is to understand how the structure of nonlinear feedbacks determines qual- 
itative behavior of the system, in particular existence of stable spatially heterogenous 
patterns. We discuss different mechanisms of pattern formation, i.e., diffusion-driven 
instability and hysteresis-driven pattern formation. The first class of models uses spe- 
cial features of diffusion, which results in the destabilization of the spatially homoge- 
neous steady state and emergence of spatial heterogeneity. The second mechanism of 
pattern formation in such systems is based on the existence of multiple steady states 
and hysteresis in the intracellular dynamics. Diffusion of the signaling molecules tries 
to average different states and is the cause of spatio-temporal patterns. 


8.1.2 Mechanisms of Developmental Pattern Formation 


One of the crucial issues in developmental biology is to understand how coordinated 
systems of positional information are established during an organism’s development 
and how cells in the organism respond to the associated signaling cues, processes 
which ultimately result in the subdivided and patterned tissues of multicellular organ- 
isms [35, 36,46]. Experiments suggest that during development cells respond to local 
positional cues that are dynamically regulated. The hypothesis is that cells differenti- 
ate according to positional information [46]. The question is how this information is 
supplied to the cells. 

There exists a number of models for pattern formation and regulation based on the 
idea that positional information is supplied to cells by a diffusing biochemical mor- 
phogen [7,37, 46]. It links the expression of target genes with local concentrations of 
morphogen molecules (ligands). Different concentrations of morphogens are able to 
activate transcription of distinct target genes and thus cell differentiation. However, 
both regulatory and signaling molecules (ligands) act by binding and activating re- 
ceptor molecules which are located in the cell membrane (or, with lipophilic ligands, 
in the cytoplasm) [20, 36]. This observation leads to a hypothesis that the positional 
value of the cell may be determined by the density of bound receptors which do not 
diffuse [35]. 
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8.1.3 Motivating Application: Pattern Control in Hydra 


One of the most frequently discussed organisms in theoretical papers on biological 
patterns formation is the fresh-water polyp Hydra. What is peculiar about Hydra? 

Hydra, a fresh-water polyp, is one of the oldest and simplest multi-cellular organ- 
isms equipped with typically animal cells such as sensory cells, nerve cells and muscle 
cells. The animal has an almost unlimited life span and regeneration capacity. Similar 
to plants, in Hydra tissue there are stem cells that are constantly dividing and regen- 
erating the adult structures of the polyps. This unlimited growth indicates that Hydra 
does not undergo senescence and, in this sense, itis biologically immortal [4]. Morpho- 
genetic mechanisms active in adult polyps are responsible for the regenerative ability 
and the establishment of a new body axis. Research on Hydra might reveal how to se- 
lectively reactivate the genes and proteins to regenerate human tissues. The fact that no 
tumor formation or other malignancies have been reported for Hydra so far, indicates 
that growth control and tissue homeostasis in normal Hydra polyps are very efficient. 

The developmental processes governing formation of the Hydra body plan and its 
regeneration are well understood at the tissue level [35,36]. Therefore, experiments 
performed on Hydra provide a good ground for testing the abilities and limitations of 
mathematical models. We may distinguish three main experiments: 


e De novo pattern formation. 
It was shown that normal Hydra can regenerate from random cellular aggregates [10, 
36] (see Figure 8.1). Reorganization does not result from a spatial rearrangement, 
but it is an effect of concerted changes in the functional state of the cells. The cells 
do not sort with respect to the positional origin along the body axis [39]. These 
experiments suggest that there exist mechanisms which define new centers of head 
organizing activity within an initially chaotic mass of cells. 


e Cutting experiments. 

Hydra has a high capacity to regenerate any lost body part, which occurs mainly 
by the re-patterning of existing tissues and is an example of morphallaxis [46]. The 
lack of growth requirements for regeneration is shown in heavily irradiated polyps. 
No cell divisions occur, but the animals can still regenerate normally. Consequently, 
the mechanism of pattern formation in Hydra seems to be independent of growth. 
Overlapping cut levels show that the same cells can form either the gastric region, 
or the head, or the foot, according to their position along the body axis (see Fig- 
ure 8.2). The experiment shows that after a transverse cut both parts of the animal 
can regenerate [35]. Moreover, the polarity is maintained even in small pieces of 
the body. A tissue piece containing 150-300 epithelial cells, i.e., about 1 percent 
of normal polyp, regenerates a complete Hydra [41]. However, below this size no 
regeneration takes place. There are also observations showing that the time required 
for the regeneration decreases with increasing tissue size [41]. 
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Figure 8.1. Time evolution of aggregates from randomly arranged Hydra cells. Several polyps 
were disintegrated into a suspension of isolated cells, which subsequently were allowed to 
reestablish contact. Within two weeks the aggregate reintegrated itself into intact animals 
[Courtesy of W. Miiller]. 
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Figure 8.2. Cutting experiment. Hydra regenerates after a transverse cut of cells of the gastric 
region (from both upper and lower half of the body column). In one experiment (left-hand 
side) the lower body column is removed, in the second experiment (right-hand side) the upper 
part is removed. The cut levels are not identical but somewhat different to show that one and 
the same group of cells (marked in grey) can form a foot (left-hand side), or a head (right-hand 
side) or a gastric segment (original state in the middle). The function of the cells depends on 
the position along the body column [Courtesy of W. Miiller]. 


e Grafting experiments. 
Grafting experiments show how disparities between the positional value of the trans- 
plant and the surrounding host tissue result in the head or foot formation leading to 
development of new organisms with multiple heads or feet [35,36] (see Figure 8.3). 
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Figure 8.3. Grafting experiment. Determination of relative positional information values by 
transplantation. Pieces of tissue are grafted from one animal to another and one of three out- 
comes is observed. (1) If the tissue is transplanted from the upper position along the body 
column to the lower position then a new head is formed. (2) If the former and new position is 
the same then the piece is integrated and nothing is observed. (3) If the tissue is grafted to the 
upper position a new foot is formed [Courtesy of W. Müller]. 
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Figure 8.4. The illustration of the idea of “positional value”, which is supplied to the cells and 
interpreted by them. The hypothesis is that the formation of the head is determined by the high 
“positional value” (which is above some threshold). The figure shows the “positional value” 
for a supernumerary head structure [Courtesy of W. Miiller]. 
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To conclude, experiments of this kind suggest that the cells respond to local positional 
cues that are dynamically regulated. It leads to the hypothesis of Wolpert on posi- 
tional information [46] (compare Figure 8.4). The question is how this information 
is supplied to the cells and which mechanisms control the formation of spatially het- 
erogenous structures in the positional information, and consequently patterns of cell 
differentiation. 

In the remainder of this paper we will address this question based on the results 
of mathematical models employing different hypotheses. Since the exact molecular 
mechanism of pattern formation in Hydra is unknown, the proposed models are hypo- 
thetical. They attempt to answer the following questions: 
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e What minimal processes are sufficient to produce de novo patterns? 


e Which models are able to capture the results of the above experiments? 


The problem was first approached by Wolpert [46], who suggested a gradient model 
to account for head formation, in which at the head end a morphogen S is emitted. 
The morphogen spreads by diffusion and is distributed down the body. This diffus- 
ing chemical induces formation of the head. The proposed model corresponds to the 
assumption that morphogens are secreted only by a group of cells in some restricted 
region of a tissue and then transported in an adjacent tissue. While there is experimen- 
tal evidence of such signaling in other systems, such as, for example, Drosophila wing 
imaginal disc or Spemann organizer, it is not the case in Hydra de novo development 
from the dissociated cells [36]. 


8.1.4 Diffusive Morphogens and Turing Patterns 


The question of de novo pattern formation in a homogenous tissue was addressed by 
Turing in his pioneering paper [45]. He proposed a hypothesis that can be stated as fol- 
lows: When two chemical species with different diffusion rates react with each other, 
the spatially homogeneous state may become unstable, thereby leading to a nontrivial 
spatial structure. 

The idea looks counterintuitive, since diffusion is expected to lead to the uniform 
distribution of the particles. Mathematical analysis of reaction-diffusion equations pro- 
vides an explanation for the phenomenon postulated by Turing. The proposed mech- 
anism of pattern formation is related to a local behavior of solutions of a reaction- 
diffusion system in the neighborhood of the constant solution that is destabilized via 
diffusion. Patterns arise through a bifurcation, which we call diffusion-driven instabil- 
ity (DDI). They can be spatially monotone corresponding to the gradients in positional 
information or spatially periodic. 


Definition 8.1 (Turing instability). A system of reaction-diffusion equations exhibits 
DDI (Turing instability) if and only if there exists a constant stationary solution, which 
is stable to spatially homogenous perturbations, but unstable to spatially heterogenous 
perturbations. 


The original idea was presented by Turing on the example of two linear reaction- 
diffusion equations of the form 


oe =DAu+Au inQ, 
ot 
d,u(t,0) = 0 on dQ, 
u(0, x) = u(x), (8.1) 


where u € RÈ is a vector of two variables, D is a diagonal matrix with nonnega- 
tive coefficients du, dy on the diagonal, the symbol 0, denotes the normal derivative 
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(no-flux condition), and Q is a bounded region. Here the only constant steady state is 
(0, 0). 
Following Turing [45], we can formulate the following result on DDI: 


Theorem 8.2 (Allan Turing). Assume that trA < 0, detA > 0 and dy > 0. There 
exists dy > 0 (small enough) such that the constant steady state (0,0) is unstable for 
the reaction-diffusion equation (8.1). 


It can be proven using a spectral decomposition of the Laplace operator with ho- 
mogenous Neumann boundary conditions and calculating the eigenvalues of obtained 
finite dimensional operator. 

Due to the local character of Turing instability, the notion has been extended in a 
natural way to the nonlinear equations using linearization around a constant positive 
steady state. However, in case of nonlinear systems we may deal with the existence of 
multiple constant steady states. In such cases we observe the existence of heterogenous 
structures far from the equilibrium and the global behavior of the solutions cannot 
be predicted by the properties of the linearized system, e.g., [27,43]. In fact we can 
observe a variety of possible dynamics depending on the type of nonlinearities. On 
the other hand, Turing instability can also be exhibited in degenerated systems such as 
reaction-diffusion-ODE models or integro-differential equations, for example, shadow 
systems obtained through reduction of the reaction-diffusion model [26, 27]. 

Following all these observations and original character of Turing’s system we define 
the Turing patterns in the following way: 


Definition 8.3. By Turing patterns we refer to the solutions of reaction-diffusion equa- 
tions that are 


e stable, 

e stationary, 

e continuous, 

e spatially heterogenous and 

e arise due to the Turing instability (DDI) of a constant steady state. 


It can happen in a reaction-diffusion system with DDI that all nonconstant stationary 
solutions are unstable and then the solution converges to another constant solution 
or to a dynamical structure such as a spike pattern [26, 43]. In case of at least three 
equations, the system can also exhibit a Turing-type Hopf bifurcation, which leads to 
spatio-temporal oscillations [19]. 
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8.1.4.1 Activator-inhibitor Model 


The most famous realization of Turing’s idea in a mathematical model of biological 
pattern formation is the activator-inhibitor model proposed by Gierer and Meinhardt 
[7]. The model aims to explain head formation in Hydra due to a coupling of a local 
activation to a long-range inhibition process. An activator promotes head formation 
and increases itself autocatalytically. An inhibitor acts as a suppressant against the self- 
enhancing activator to prevent the system from unlimited growth. In this approach the 
positional value is interpreted as the density of the activator. Gradients of morphogens 
are formed by the DDI mechanism. Each of the various body parts is assumed to be 
under control of a separate activator-inhibitor system (for details, see [32]). The basic 
activator-inhibitor model takes the form 


ð 3? a? 

le = Da age T Pa + Oa — Haa, 

ð 3? 5 

y = Dryl + Pha? + Oh — Hph, (8.2) 


where a and h denote the concentrations of the activator and the inhibitor, respectively. 
The parameters og and op describe de novo production, Ha and up are the rates of 
degradation and pg and pp the parameters of the activator-inhibitor interactions. The 
model and several of its modifications were applied in the study of various topics from 
developmental biology (see, e.g., [33,34]). Due to its interesting mathematical features 
and emerging singularities, the model has also attracted a lot of attention from the side 
of mathematical analysis, e.g., [31,43]. 

The activator-inhibitor theory operates with purely hypothetical morphogens. As we 
can see in Theorem 8.2, the key mechanism of Turing-type patterns is that an inhibitor 
diffuses faster than an activator. However, dynamics and complex tissue topologies 
are likely to prevent the establishment of long-range inhibitor gradients. Furthermore, 
diffusion rates of typical morphogens are often found to be quite small [9], i.e., do not 
allow significantly varying diffusion rates as required by the Turing mechanism. In 
the case of Hydra, while recently Wnt can be identified as an activator [10], a long- 
range inhibitor is missing [21,34]. These observations support the search for a dif- 
ferent inhibitory mechanism such as mechanical inhibition [6] or different than DDI 
mechanism of pattern formation [24]. In the context of Hydra experiments it is also 
important to note that the shape of Turing patterns depends on the size of the domain 
and diffusion rather than on initial conditions. Therefore, one of the difficulties of the 
Turing-type models is their inability to reproduce the experiments resulting in multiple 
head formation in Hydra [23]. 
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8.1.5 Receptor-based Models 


Another type of mathematical models for pattern formation follows the hypothesis that 
the positional value of the cell is determined by the density of cell-surface receptors, 
which regulate the expression of genes responsible for cell differentiation [35], see 
Figure 8.5. 
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Figure 8.5. Bound receptors density determining “positional value”: the head is formed if the 
density of bound receptors is high (above some threshold). Consequently, in normal develop- 
ment we expect a gradient-like distribution of bound receptors [Courtesy of W. Müller]. 


The receptor-based models are based on the idea that epithelial cells secrete ligands 
(a regulatory biochemical), which diffuse locally within the interstitial space and bind 
to free receptors on the cell surface [23,29]. It results in a bound receptor that can be 
removed from the cell surface due to degradation or internalization, or dissociate back 
to free receptors and ligands. Both ligands and free receptors are produced within the 
whole tissue and undergo natural decay. 

The first receptor-based model for Hydra was proposed by Sherrat, Maini, Jäger 
and Müller in [40] in the following form (SMJM model), 


ð 3? 
ap? = Pega. + Sq(x) — Haa — keae — kaaf + kab, 


af = kab —kaaf + kilo(x) + Bb — f1, 


ð 
ape T Kaaf — (ka + ki)b, 


ə 3? 
E = Dega? + Se(x) — Hee, (8.3) 


defined on a bounded one-dimensional domain with zero-flux boundary conditions for 
a and e. The variables f, b, a and e denote the density of free receptors, bound recep- 
tors, biochemical (ligands) and enzyme, respectively. In order to achieve the required 
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local competition phenomenon, it is assumed in the SMJM model that the terms de- 
scribing de novo production of free receptors, ligands and enzyme depend on the po- 
sition of the tissue along the body axes x € [0, L] and also on y(x), which denotes 
position from which the tissue at location x originates. The functions describing the 
production of new receptors and production of enzyme are linearly decreasing in x. 
The production of ligands is assumed to be constant on 4/5 of the domain and then 
decrease linearly to zero. It is assumed that 


a(x) = ay[1 — y(x)/L] + æ2y(x)/L, 
Se(x) = si[1 — y(x)/L] + soy(x)/L, 


where sq(x) is constant for y(x) € [0, 44] and decreases linearly to zero for y(x) € 


[(S, 1]. The combination of these two parallel gradients enables the model to cap- 
ture some results of grafting and cutting experiments. Thus, the model functions not 
because of nonlinear interactions between receptors and ligands, but because of the 
assumption that cells produce new molecules depending on the position they had in 
the donor organism. In conclusion, the SMJM model is not a model for de novo pattern 
formation. 

Later, receptor-based models without imposing initial gradients were proposed by 
Marciniak-Czochra [23]. In general, equations of such models can be represented by 
the following initial boundary-value problem, 


us = DAv+ f(u, v) ing, 


vs = g(u, v) in Q, 
dnu = 0 on 0Q, 
u(x, 0) = uo(x), v(x, 0) = vo(x), (8.4) 


where u is a vector of variables describing the dynamics of diffusing extracellular 
molecules and enzymes, which provide cell-to-cell communication, while v is a vec- 
tor of variables localized on cells, describing cell surface receptors and intracellular 
signaling molecules, transcription factors, mRNA, etc. D is a diagonal matrix with 
positive coefficients on the diagonal, the symbol 0, denotes the normal derivative 
(no-flux condition), and Q is a bounded region. 

A rigorous derivation, using methods of asymptotic analysis (homogenization) of 
the macroscopic reaction-diffusion models describing the interplay between the non- 
homogeneous cellular dynamics and the signaling molecules diffusing in the intercel- 
lular space has been undertaken in [25, 29]. It is shown that receptor-ligand binding 
processes can be modeled by reaction-diffusion equations coupled with ordinary dif- 
ferential equations in the case when all membrane processes are homogeneous within 
the membrane. 

As shown in [23] and also more recently highlighted in [16], receptor-based mod- 
els may exhibit Turing-type instability. The simplest receptor-based model takes into 
account only one type of diffusive signaling molecules. The basic model of this type, 
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for one-dimensional epithelial sheet, takes the form 


0 
crt ie = —HUfrf + pr(rf,rp)—brfl + drp, 


ð 
a? = —pprp + brfl — drp, 

ae ae ee 1+ pi(l.rp) +d (8.5) 
=l = — — br r rb, : 
at y ox 7 Hı f Pl b b 


with zero flux boundary conditions for l, 0,/(t,0) = 0;/(t, 1) = 0. 

The model takes into account dynamics of free and bound receptors on cell 
membranes, denoted by r¢(t,x) and rp(t, x), respectively, and diffusing signaling 
molecules denoted by /(t, x). The original model operated with purely hypothetical 
molecules. However, in case of Hydra pattern formation we may associate the vari- 
ables to Frizzled receptors and Wnt ligands [17]. The hypothesis that Wnt is a diffusing 
ligand is supported by experimental evidence [10, 21]. 

The effects of intracellular dynamics are modeled via nonlinear functions describing 
production of new signaling molecules and free receptors, p; and p,, respectively. 
Besides, the kinetics describe binding at a rate b, dissociation at the rate d and natural 
decay at the rates uf, pp and py. y = Z is a scaling coefficient depending on the 
domain length L and the diffusion oa at dj. 

As it was stated in [28, Proposition 3.1], a generic system of two ordinary differential 
equations coupled with a reaction-diffusion equation exhibits DDI if there exists a 
positive, spatially constant steady state, for which the following conditions are satisfied 


—tr (A) > 0, (8.6) 

—tr (A) È` det (Ajj) + det (A) > 0, (8.7) 
i<j 

—det (A) > 0, (8.8) 

—det (A12) > 0, (8.9) 


where A = (41); j=1,2,3 
earized around this constant positive equilibrium and A;; is a submatrix of A consisting 
of the 7 th and j th column and i th and j th row. 

Conditions (8.6)—(8.8) are necessary for the stability of the steady state in the ab- 
sence of diffusion. Inequality (8.9) is a sufficient and necessary condition for destabi- 
lization of this steady state. These conditions guarantee that the model (8.5) exhibits 
diffusion-driven Apay if the function p, evaluated at the steady state satisfies 
Prf, rb) < Ffa a Pr(Tf, Tp), which, in particular, holds for p, (rf) = mırt! for 


is the Jacobian matrix of the system without diffusion lin- 


a > 0. This meguality can be interpreted as an autocatalysis at the steady state in the 
first equation of the system. This condition leads to the instability of those constant 
solutions, for which an autocatalysis occurs. 
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Figure 8.6. Spatial profile of rp for different initial perturbations and a fixed y. On the left- 
hand side we present the initial condition and on the right-hand side the final pattern origi- 
nating from such an initial condition. The perturbation in initial data is of order 107? (it can 
be arbitrarily small) while the final peak is of height 8 x 10*. (a) Different initial conditions 
and corresponding solutions are depicted using matching line styles. The location of the peak 
strongly depends on the initial condition. (b) The results for the initial conditions with two 
maxima or two minima—the result is always one peak (depicted using matching line styles). 


Following the classical Turing idea, one expects stable patterns to appear around 
the constant steady state in the system with DDI property. Interestingly, in numeri- 
cal simulations of model (8.5), diffusion-driven instability of the constant steady state 
leads to the emergence of growth patterns concentrated around discrete points along 
the spatial coordinate, which take the mathematical form of spike-type spatially inho- 
mogeneous solutions [23]. The structures are not robust and depend strongly on initial 
conditions. In some cases, blow up occurs. Definitely, the observed solutions are not 
Turing patterns, see Figure 8.6. 

Recent analytical studies of the reaction-diffusion-ODE models with only one 
nonzero diffusion coefficient revealed that monotone or periodic stationary solutions 
can be constructed for most interesting models [8, 26,27]. However, the same mecha- 
nism that destabilizes constant solutions of these models also destabilizes non-constant 
solutions [26]. Consequently, there exist no stable continuous stationary solutions for 
the initial boundary-value problems for ordinary-PDE systems as the one in (8.5). 
While in some other applications such dynamical spike patterns are of biological rel- 
evance [27, 28, 38], they cannot be applied to describe pattern formation in Hydra. 

Considering two diffusing signaling factors leads to a four-variable receptor-based 
model exhibiting Turing patterns [23]. In this model it is additionally assumed that 
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there exists a second diffusing substance, functioning as an enzyme as in the SMJM 
model [40], which is secreted by cells, diffuses along the body column and degrades 
the ligands. The equations have the following form 


0 
alr = er (rp) + prls: te) — (rp. 1) + dr), 


ð 
Tb = —Hb (rp) + b(rf, 1) — d(rp), 


ðt 
0 02 
ar! = dina! — wi) — b(rf, l) + pire. ro) + dro) — belle), 
0 32 
a = dex ae — Mele) + Pell. ro), (8.10) 


with zero flux boundary conditions for / and e. Here e denotes the density of enzyme, 
be the rate of binding of ligands and enzyme, pe the rate of production of enzyme, 
Le the rate of decay of the enzyme, de is the diffusion coefficient for enzyme, and the 
other terms are as in the three-variable model. 

The role of the enzyme is to remove the biological regulator, i.e., ligand, before it 
binds to the receptors on the cell surface. It is important to stress that this model can- 
not be simplified to an activator-inhibitor system of the type (8.2). The four-variable 
receptor-based model consists of two subsystems. The reaction-diffusion subsystem 
describing ligand and enzyme dynamics is not of the activator-inhibitor type and can- 
not produce diffusion-driven patterns itself. It is the ODEs subsystem that causes desta- 
bilization of a constant solution and emergence of a Turing pattern. 

In such a model, patterns can evolve due to the DDI, even if no self-enhancement 
of free receptors nor ligands is assumed. p, is assumed to be a function of rp, since 
it is known from a number of other biological contexts that there can exist a posi- 
tive feedback loop between the density of bound receptors on the cell surface and the 
subsequent expression of new receptors [15,44]. Also no assumption on the range of 
enzyme diffusion is needed. 

These observations show that including receptor dynamics in the model of inter- 
acting diffusing signaling molecules allows to relax the assumptions on the range of 
diffusion and type of nonlinear interactions necessary for a formation of stable spa- 
tially heterogenous patterns. In particular, it seems that although the Wnt antagonist 
Dickkopf found in Hydra tissue [2] does not satisfy the assumptions of the inhibitor 
from the Gierer-Meinhardt model [34], its interactions with Wnt signaling may lead 
to a stable gradient-like pattern formation as in the receptor-based model. 

In the four-variable receptor-based model, similarly to the activator-inhibitor model, 
the spatially homogeneous steady state bifurcates into the spatially inhomogeneous 
solution which has a maximum at the one end and a minimum at the other end for some 
range of the domain size. This model is robust and the final pattern does not depend on 
small perturbations of initial conditions. Pattern formation phenomenon is similar to 
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Figure 8.7. Numerical simulation of the cutting experimentin the four-variable receptor-based 
model. Left-hand panel: Initial data (dashed line) corresponding to a surgical removal of the 
lower part (half) of the body column. We observe that a reorganization of the “gradient” on a 
smaller domain corresponds to the formation of a new “foot”. Right-hand panel: Initial data 
(dashed line) corresponding to a surgical removal of the upper part (half) of the body column. 
We observe that a reorganization of the “gradient” on a smaller domain corresponds to the 
formation of a new “head”. 


that in the activator-inhibitor model of Gierer and Meinhardt (8.2) with the difference 
that maxima of the pattern have the shape of plateaux and not spikes as was the case in 
(8.2). It is related to uniform boundedness of solutions in model (8.10). Models with 
Turing patterns can describe self-organization of Hydra cells (see Figure 8.8, left-hand 
panel) and are able to simulate the cutting experiments (see Figure 8.7). Concluding, 
introduction of a second diffusing biochemical species improved performance of the 
model. The four-variable receptor-based model can explain at least as much as the 
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Figure 8.8. Left-hand panel: De novo gradient-like pattern formation in four-variable receptor- 
based model. Right-hand panel: Simulation of a transplantation experiment for the initial data 
corresponding to the head grafting. The final distribution shows the transplant disappearance. 
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activator-inhibitor model regarding de novo pattern formation and basic experiments. 
Numerical studies of both models based on Turing mechanism showed that the grafting 
experiments could not be explained within such an approach without changing the size 
of domain (or diffusion coefficient), which does not reflect experimental conditions 
(see Figure 8.8, right-hand panel). 


8.1.6 Multistability 


Transplantation experiments suggest that there are a number of locally asymptotically 
stable patterns depending on the past history. The patterns may have multiple peaks, 
which depend on the local cues induced by the grafted tissue. Such experiments sug- 
gest a mechanism of pattern formation based on multistability in intracellular signal- 
ing. Coupling diffusion with a kinetics system with multiple stable steady states and 
hysteresis may lead to the coexistence of different patterns for the same parameters 
but depending on the initial conditions. Such a hypothesis was incorporated in the 
receptor-based model proposed in [24] by replacing the function p;, describing the 
rate of production of diffusing signaling molecules, by a new variable modeled using 
an additional ordinary differential equations. The model includes a hysteresis-based 
relation in the quasi-stationary state in the ODEs subsystem, i.e., g(u, v) = 0 in the 
system of equations (8.4) (see Figure 8.9). 


Figure 8.9. A typical configuration of the kinetic functions in a receptor-based model (8.4) 
with hysteresis in the quasi-stationary ODEs subsystem. 


The model suggests how the nonlinearities of intracellular signaling may result in 
spatial patterning. It allows for formation of gradient-like patterns corresponding to the 
normal development as well as emergence of patterns with multiple maxima describ- 
ing transplantation experiments (see Figure 8.10). Numerical simulations show the 
existence of stationary patterns resulting from the existence of multiple steady states 
and switches in the production rates of diffusing molecules, see Figure 8.10, left-hand 
panel. The patterns observed in such models are not Turing patterns. In fact, the system 
does not need to exhibit DDI. Indeed, in most cases its constant steady states do not 


Bereitgestellt von | De Gruyter / TCS 
Angemeldet 
Heruntergeladen am | 16.10.19 13:17 


8.1 Reaction-Diffusion Models of Pattern Formation in Developmental Biology 207 


80+ 
gy 807 p 
2 6014 = a 
2 Foot 2 
® 40 2 
o oT 
3 0 > 3 
3 10008 


Ox 
0 


00 Time 


100 9 


Figure 8.10. Simulations of the receptor-based model with hysteresis. Formation of a gradient- 
like pattern corresponding to a normal development and head formation in Hydra (left-hand 
panel) and formation of two heads pattern for the initial conditions corresponding the trans- 
plantation experiment (right-hand panel). 


change stability. In such models, spatially heterogenous stationary solutions appear 
far from equilibrium due to the existence of multiple quasi-steady states. 

Properties of a hysteresis-based mechanism of pattern formation have been recently 
studied in a minimal version of the model consisting of one reaction-diffusion equation 
coupled to one ODE [17]. In such a model, infinitely many stationary solutions can be 
constructed. Such solutions are discontinuous in the nondiffusing variables. The shape 
of spatial structures can be very irregular, since it depends strongly on the initial con- 
ditions. Therefore, the model can simulate the effects of transplantation experiments 
and formation of multiple heads, see Figure 8.10, right-hand panel. 

On the other hand, it was shown that the system with multistability but reversible 
quasi-steady states in the ODE subsystem, i.e., g (u, v) globally invertible, cannot ex- 
hibit stable spatially heterogeneous patterns. Hysteresis is necessary to obtain stable 
patterns. 


8.1.7 Discussion 


Transplantation and tissue manipulation experiments provided data for models of pat- 
terning in Hydra, starting with the positional information ideas of Wolpert [46], the 
activator-inhibitor model of Gierer and Meinhardt [7, 32-34] and, finally, receptor- 
based models of Marciniak-Czochra [17,23,24]. Each model has shed light on differ- 
ent but overlapping aspects of self-organization and regeneration. Now, it is possible 
to state which conceptual elements have to be present in a complete model, although 
modeling of Hydra regeneration still involves quite a few unsolved problems. In the 
framework of reaction-diffusion systems there are essentially two ways in which a 
system of identical cells can start to differentiate: 
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e There is a critical number of cells (size of domain), above which the spatially homo- 
geneous attractor loses stability, which leads to “spontaneous” spatial patterning. It 
is the case for the models with the Turing instability. Such models can explain de 
novo pattern formation since for some set of parameters and the domain size value, 
the final pattern is the same and does not depend on the initial perturbation. 


e There is an external inducing signal which drives the system into a new, spatially 
inhomogeneous state. Such a signal originates from another group of already differ- 
entiated cells and it must be strong enough to trigger differentiation. It corresponds 
to a sufficiently strong initial perturbation of the homogeneous steady state. This 
type of initialization of the pattern-forming mechanism is involved in the model 
with hysteresis. 


The experiments showing de novo formation of Hydra from the dissociated cells 
could suggest a Turing-type mechanism. On the other hand, transplantation exper- 
iments suggest coexistence of different spatially inhomogeneous stationary patterns 
which grow up for different initial conditions. Experiments show that large pertur- 
bations (but within the range of the values of the solution itself) of the gradient-like 
solution should lead to another solution. The observations may be explained by mod- 
els with multiple steady states exhibiting hysteresis. In such models solutions depend 
on the initial condition similarly to the Hydra resulting from the grafting experiment 
depends on the graft position and not on the size of the animal. 

The question is whether these two kinds of experiments can be explained using the 
same mechanism and whether it could be the combination of the already considered 
mechanisms. To clarify these issues new models including a more detailed description 
of cell-to-cell and intracellular signaling should be developed. Mathematical under- 
standing of the relation between the structure of nonlinearities and the dynamics of 
model solutions shall be helpful both in building new models and also in designing 
experiments that might help to verify different hypotheses. 

Many uncertainties exist regarding the biological foundations of the models. Fur- 
ther biological discoveries are needed to gain insight into the molecular nature of cell 
communication and of the positional value. To understand the morphogenesis in Hy- 
dra it is necessary to bridge the gap between experimental observations at the cellular 
level and those at the genetic and biochemical levels. 

The advent of new techniques in molecular biology has recently made it possible 
to advance the understanding of the development of multicellular organisms. Large 
scale expression screening helps to identify new factors involved in embryonic devel- 
opment. Recently, expression analysis during regeneration and budding indicated a 
pivotal role of the Wnt (wingless gene) pathway in the Hydra head organizer [10]. Also 
the evidence of Dkk (Dickkopf) signaling in Hydra regeneration was provided [2]. Ex- 
perimentally observed patterns of Wnt and Dkk gene expression give rise to many new 
questions. New aspects of Wnt-Dkk signaling, such as bi-stability in Wnt dynamics 
and switches in the Dkk functionality depending on the cellular context, were also 
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recently found in other model organisms [14, 18,22]. Mathematical modeling should 
integrate these concepts and observations into a new model of pattern formation con- 
trolled by the intracellular dynamics of Wnt-Dkk signaling. 

In conclusion, growth and pattern formation provide a great source of interesting 
and novel mathematical problems, while mathematics can be used as a tool to explore 
different mechanisms and processes underlying these phenomena. The use of realistic 
models may help to understand many complex processes. 
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